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I.  INTRODUCTION  AND  SUMMARY 


In  support  of  POET  analysis  of  the  detection  and  acquisition  ranges  of  a  two-color 
sensor  on  a  terminal-based  interceptor,  a  procedure  has  been  set  up  to  calculate  the  transient 
temperature  distribution  on  the  surface  of  a  reentry  vehicle  (RV).  The  example  exercise  of 
the  procedure  provided  here  is  for  an  RV  with  a  reference  ballistic  coefficient  of  about  1500 
lb/ft2  reentering  at  a  reentry  angle  near  25  deg  on  a  ballistic  path  with  a  range  of  about 
10,000  km  to  impact  (a  reentry  velocity  of  about  23,150  ft/sec).  The  inputs  to  the  IDA 
RANGE  trajectory  program  that  were  used  to  characterize  the  hypothetical  RV  are  the 
following: 

Cone  half  angle  =  10  deg. 

Ratio  of  hemisphere  cross-sectional  area 
to  reference  area  =  0.041, 

Ratio  of  base  area  to  reference  area  =  1.00, 

Ratio  of  surface  area  to  reference  area  =  5.76, 

Mass  =  750  lb,  and 

Reference  area  =  4.55  ft2. 

These  numbers  are  consistent  with  a  nose  radius  of  about  3  in.  (0.077  m)  and  an 
overall  length  of  about  5.9  ft  (1.8  m).  The  reference  ballistic  coefficient  is  in  effect  defined 
using  the  hypersonic  inviscid  drag  coefficient,  i.e.,  that  for  velocities  above  about  Mach  10 
where  the  base  drag  coefficient  has  diminished  to  negligible  values  and  for  altitudes  below 
about  200  kft  where  the  skin-friction  drag  coefficient  has  diminished  to  negligible  values 
(IDA,  1969).  At  lower  Mach  numbers  and  at  higher  altitudes,  the  overall  drag  coefficient 
generally  increases  above,  and  the  ballistic  coefficient  decreases  below,  the  reference  value. 
The  reentry  angle  is  defined  as  the  path  angle  measured  from  the  horizontal  at  an  altitude  of 
400  kft. 

The  results  of  the  example  calculations,  considering  free-molecule  or  continuum- 
flow  (laminar)  aerodynamic  heating  in  the  appropriate  altitude  regimes  (with  "bridging"  in 
the  transition  between  them),  reradiation  with  an  emissivity  of  0.25,  and  one-dimensional 
transient  heat  conduction  into  a  glass-fiber-phenolic  heat  shield,  indicate  that  the  stagnation 
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point  on  the  nose  teaches  the  "ablation"  temperature,  arbitrarily  taken  as  2000  K,  at  an 
altitude  of  about  260  kft,  and  that  the  surface  of  all  of  the  conical  body  passes  through  a 
temperature  of  1000  K  at  an  altitude  of  about  250  kft  and  reaches  the  "ablation"  temperature 
at  an  altitude  of  about  200  kft 
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II.  ANALYSIS 

The  principal  questions  addressed  by  the  analysis  were  the  following: 

(1)  What  governs  the  input  heating  at  different  points  on  the  RV? 

(2)  What  governs  the  thermal  response  of  the  RV  surface? 

Considerations  involved  in  (1)  include  the  determination  of  (a)  the  functional 
dependences,  and  altitude  regions  of  applicability,  of  ffee-molecule  heating  and  laminar 
convection  at  the  stagnation  point,  (b)  the  flight  conditions  on  the  particular  reentry 
trajectory  as  a  function  of  time,  and  (c)  the  heating  distribution  at  points  other  than  the 
stagnation  point  on  a  hemisphere-cone  body. 

Considerations  involved  in  (2)  include  determinations  of  (a)  the  transient  thermal 
response,  i.e.,  time-dependent  heat  conduction  into  the  heat-shield  material,  (b)  typical 
physical  properties  of  a  glass-fiber-phenolic  ablative  heat  shield,  and  (c)  the  dependence  of 
surface  temperature  on  emissivity  (a  discretionary,  not  inherent,  property  of  the  surface). 

The  heat  transfer  rate  at  the  stagnation  point  in  the  free-molecule  flow  regime,  qpM, 
i'  given  in  Gilbert  and  Scala,  1965,  as 


wfp-vi  i  + 


im  s„ 


ft-lb/(ff  sec) 


where 


a  =  accommodation  coefficient  (taken  as  1.0) 

poo  =  local  free-stream  atmospheric  density  (slug/ft3) 

Voo  -  vehicle  flight  velocity  with  respect  to  the  local  atmosphere  (ft/sec) 

Sop  =  Mjy/2 

Moo  =  vehicle  flight  Mach  number 

Y  =  ratio  of  the  specific  heats  of  air. 
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For  Mach  numbers  of  16  or  greater  (typical  of  the  heating  regime  of  interest  for  reentry- 

vehicle  detection  and  tracking),  the  molecular-speed  correction  term  [  \}(ITk  sj)  has  a 
value  less  than  two  percent 


Ignoring  the  molecular-speed  correction  term,  and  changing  the  Units  to  more 
convenient  (for  comparison  with  laminar  heat-transfer  below)  values,  the  formula  for  free- 
molecule  heating  at  the  stagnation  point  becomes 


qFM  =  1.528  x  10* 


(  Po.  ) 

lp°  J 

JO4, 

BtuZ(frsec)  , 


where  po  is  the  atmospheric  density  at  sea  level. 


The  here-adopted  empirical  correlation*  for  the  hypersonic  laminar  ("continuum") 
heat-transfer  rate  at  the  stagnation  point,  qL,  is  (per  IDA,  1966,  p.  22,  leaving  out  for  the 
moment  a  factor  depending  on  surface  temperature  which  is  close  to  1  for  most  cases  of 
interest) 


865 


,3.15 


J  Btu/(ft2  sec) 


where  Rn  is  the  hemispherical-nose  radius  in  feet.  The  missing  correction  factor,  involving 
the  flight  velocity  and  a  value  of  the  yet-to-be-determined  surface  ("wall")  temperature  Tw, 
is 

vj-  6.75  x  106  CIV289) 

\d-6.75xlO* 

where  Tw  is  in  degrees  Kelvin.  This  factor  will  be  included  later  in  the  transient  analysis, 
in  which  a  value  of  the  wall  temperature  becomes  available. 

The  above  two  formulas  for  the  stagnation-point  heating  in  the  free-molecule  regime 
(qFM>  varying  as  p^)  and  the  laminar  flow  regime  (qL,  varying  as  JpZ)  yield  the  same 
value  at  a  "crossover"  altitude,  defined  by  a  "crossover"  atmospheric  density  pc  given  by 

p</p„  =(2.0230  x  10‘S/R,)vi'5  • 


This  correlation  is  in  numerical  agreement  with  that  of  Detra,  Kemp,  and  Riddell  as  validated  with  data 
in  Perini,  1975. 
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The  crossover  altitude  for  an  Rn  of  3  in.  and  a  of  23,000  ft/sec  (typical  of  the  example 

ballistic  missile  trajectory)  is  about  303  kft  for  the  1962  Model  Atmosphere  used  in 
RANGE. 

In  the  transition  regime  between  free-molecule  and  laminar  heating  the  effective 
heating  q  is  assumed  to  vary  in  a  smooth  way,  without  a  sudden  discontinuity  in  slope,  and 
is  typically  faired  between  the  two  regimes  with  a  "bridging"  function  (see  Matting,  1971). 
In  this  analysis  we  adopt  a  general  bridging  relation  that  preferentially  weights  the  lesser  of 
the  two  values  according  to 


where 

^fm^l  =  V  P«>/Pc 

A  plot  of  the  transition  behavior  of  q/qc  (where  qc  is  the  value  from  either  heat-transfer 
relation  when  equals  pc)  as  a  function  of  Poo/Pc  for  this  general  bridging  relation,  for 

values  of  n  of  1  and  2,  is  shown  in  Figure  1,  with  transition  using  the  Matting  (op.  cit.) 
bridging  function  included  for  comparison.  A  value  of  2  for  the  exponent  n  is  chosen  in 
subsequent  analysis  here;  some  other  value  (e.g.,  1.6)  could  be  easily  substituted  when 
better  data  (than  those  included  in  Matting,  1971)  become  available. 

If  the  convective  heat-transfer  rate  q  is  expressed  as  an  equivalent  radiation  rate  of  a 
surface  with  an  emissivity  e  =  1  and  no  thermal  inertia  (or  no  heat  conduction),  an 
equivalent  radiative-equilibrium  surface  temperature,  Teq,  can  be  used  to  characterize  the 
aerodynamic  heating.  The  following  formula  to  give  this  Teq  (in  degrees  Kelvin)  from  local 
flight  conditions  was  derived  from  the  above  equations  for  q  (with  n  =  2),  and  has  been 
incorporated  in  the  RANGE  trajectory  program: 


The  output  value  of  this  Teq  from  RANGE  is  used  to  describe  the  aerodynamic 
heating,  as  crTeq4  [o  =  5.672  x  10-12  watt/(cm2  K4)],  in  the  subsequent  transient  heat- 
conduction  analysis,  in  conjunction  with  the  wall-temperature  correction  factor  that  uses  the 
corresponding  output  value  of  flight  velocity  from  RANGE.  The  aerodynamic  heating,  in 
terms  of  this  "inertialess"  radiative-equilibrium  stagnation-point  temperature  (consistent 


(Velocity  Constant;  No  Wall  Temperature  Correction) 


with  an  assumed  wall  temperature  of  289  K)  for  the  desired  reentry  trajectory  is  shown  in 
Figure  2.  The  flight  velocity  for  the  wall-temperature  correction  factor  (incorporated  in  the 
transient  analysis  where  the  wall  temperature  is  determined)  varies  only  between  extremes 
of  22,960  and  23,440  ft/sec  in  the  descent  from  547  kft  to  127  kft  in  the  example  reentry 
trajectory. 

At  points  other  than  the  stagnation  point  on  a  hemispherical  nose,  the  heat  transfer 
rate  q  drops  off  with  angle  around  the  hemisphere  from  the  stagnation  point  as  shown  in 
Figure  3.  (The  influence  on  heat  transfer  at  a  downstream  point  by  the  evolution  of  gases 
from  an  upstream  point  into  the  boundary  layer  is  not  considered.)  The  heat  transfer  on  the 
conical  surface  is  assumed  to  be  constant  at  the  value  at  the  shoulder,  based  on  the 
constancy  of  pressure  on  a  conical  surface  (NAVWEPS,  1961,  p.  154)  and  the  one-to-one 
correspondence  of  heat  transfer  with  pressure  (Detra  and  Hidalgo,  1961,  Fig.  1  therein). 
From  the  curve  in  Figure  3,  the  relative  heat-transfer  rate  q/qs  (where  qs  equals  the  q 
above)  as  a  function  of  distance  s  from  the  stagnation  point  for  our  0.077-m-radius- 
hemisphere/10-deg-cone  body  is  given  in  the  following  table: 


s/R.  deg 

a/os 

s*_m 

0 

1.0 

0.000 

(stagnation  point) 

20 

0.93 

0.027 

40 

0.72 

0.054 

60 

0.45 

0.081 

80 

0.22 

0.108-1.858 

(shoulder  and  conical  surface) 

The  transient  behavior  of  the  surface  temperature  is  computed  by  the  TRIDE 
program  (IDA,  1974)  from  the  input  heating-rate  history  derived  from  the  RANGE 
inertialess  radiative-equilibrium  stagnation-point  temperatures  (Fig.  2).  The  TRIDE 
computer  program  (reproduced  in  the  Appendix)  is  a  numerical  solution  of  the  one¬ 
dimensional*  time-dependent  heat-conduction  equation 

9T_  k  dh 
dt  cp  3x2 


The  adequacy  at  one-dimensional  versus  two-dimensional  analysis  will  be  addressed  by  comparing  the 
magnitudes  of  the  gradients  in  the  two  directions,  normal  to  the  surface  and  parallel  with  the  surface.  If 
the  gradient  (i.e.,  heat  flow)  along  the  surface  is  small  in  comparison  with  the  gradient  into  the  surface, 
its  neglect  in  the  heat  flow  analysis  is  acceptable. 
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Figure  2.  Aerodynamic  Heating  Rate  at  Stagnation  Point  (Rn  =  0.077  m 
in  Terma  of  Equivalent  "Inertialeaa"  Radiative-Equilibrium  Temperature 
(for  Tw  -  289  K)  (for  ICBM  Reentry  Vehicle:  p  -  1500  ib/ft*; 
y  *  -  24.8°;  Rlmpact  *  10,020  km) 


i 


i 
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l 
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using  the  implicit-differencing  technique  (see  Appendix).  The  calculation  treats  an  infinite 
plate  of  finite  thickness  heated  on  one  side  and  insulated  on  the  other.  At  each  time  step 
(one-second  interval)  the  aerodynamic  heating,  in  terms  of  equivalent  radiation  at  the 
inertialess  temperature  Teq  corrected  with  a  surface  ("wall")  temperature  factor,  is  absorbed, 
some  is  reradiated  according  to  the  surface  temperature  Tw  and  the  emissivity  e,  and  the 
remainder  is  conducted  inward  to  raise  Tw  and  the  internal  temperature,  according  to  the 
boundary  condition 


k2L= 

dx 


-oT, 


eq 


vj-6.75x  10*0^289) 
vi-  6.75  x  106 


+  « 


:oT 1 


The  Tw  used  in  the  reradiation  term  and  in  the  correction  factor  is  extrapolated  to  the  present 
instant  by  means  of  a  cubic  extrapolation  from  the  four  previous  values.  The  temperature 
within  the  thickness  of  the  plate  is  computed  at  10  points,  the  centers  of  10  equal 
subdivisions,  e.g., 


insulated 
back  face 


The  thickness  chosen  for  the  calculations  here  is  small  enough  to  provide  adequate  detail  in 
the  temperature  contour  but  large  enough  to  include  all  the  material  that  may  experience  a 
temperature  rise  during  the  period  of  interest.  The  thermal  conductivity  k,  heat  capacity  c, 
density  p,  and  emissivity  e  are  constants  in  the  calculation.  When  the  material  in  the  first 
cell  reaches  the  "ablation"  temperature,  taken  arbitrarily  as  an  even  2000  K,  the  calculation 
is  stopped. 

The  physical  properties  for  "glass-fiber  phenolic"  for  purposes  of  this  introductory 
exercise  were  taken  from  the  Handbook  of  Materials  Science,  Vol.  HI,  p.  34,  and  are  the 
following  (temperature  dependences  were  ignored): 

k  =  0.20  Btu/(hr  ft  F)  =  0.000826  cal/(cm  sec  C), 
c  =  0.30  Btu/(lb  F)  =  0.30  cal/(gm  C),  and 
p  =  1.825  gm/cm3. 

10 
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For  these  values  and  the  experienced  heating  rates,  a  thickness  of  0.5  cm  is  found  to  be 
appropriate;  i.e.,  the  back-face  temperature  at  0.5-cm  depth  starts  to  rise  at  about  the  time 
that  the  front-face  temperature  reaches  the  assumed  temperature  of  ablation. 

The  emissivity  is  a  parameter  of  choice  by  the  RV  designer  depending  on  his 
selection  of  an  external  coating,  and  is  varied  in  the  family  of  calculations  over  the  range  of 
0.25  to  1.0. 
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III.  CALCULATIONAL  PROCEDURE 


As  a  recapitulation  of  the  analysis  described  above,  the  calculational  procedure  is 
made  up  of  the  following  steps: 

1 .  Convert  flight-vehicle  characteristics  to  RANGE  inputs. 

2.  Run  RANGE  to  provide  the  stagnation-point  aerodynamic  heating,  expressed 
as  an  equivalent  inertialess  radiative-equilibrium  temperature,  and  the  flight 
velocity  as  functions  of  time  and  altitude. 

3 .  Select  appropriate  values  of  conductivity,  heat  capacity,  density,  and  thickness 
of  the  vehicle’s  heat  shield.  (The  material  properties  can  be  temperature 
dependent.) 

4.  Insert  heat-shield  properties  and  second-by-second  values  of  the  stagnation- 
point  aerodynamic  heating,  i.e.,  the  inertialess  radiative-equilibrium 
temperature,  and  flight  velocity  in  TRIDE. 

5.  Select  values  of  the  ratio  to  the  stagnation-point  heating  for  the  heating  at 
desired  points  on  the  vehicle. 

6.  Run  TRIDE  to  give  the  time  dependence  of  the  temperature  profile  within  the 
wall  of  the  heat  shield  for  the  heat-transfer  ratios  at  the  desired  points.  Adjust 
the  thickness  of  the  heat  shield  so  that  the  temperature  of  the  back  face  would 
just  start  to  rise  at  the  end  of  the  period  of  interest,  and  return  to  step  4,  above. 
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IV.  EXAMPLE  RESULTS 


The  temperature  history  on  the  RV  surface  at  the  stagnation  point  and  on  the 
shoulder  and  conical  surface,  calculated  by  the  procedure  outlined  above,  is  plotted  as  a 
function  of  altitude  and  time  in  Figure  4.  For  a  surface  emissivity  of  0.25,  the  temperature 
at  the  stagnation  point  passes  through  1000  K  at  an  altitude  of  about  290  kft  and  reaches  the 
"ablation"  temperature  (taken  as  2000  K)  at  about  260  kft;  the  temperature  at  the  shoulder 
and  on  all  the  conical  surface  passes  through  1000  K  at  about  250  kft  and  reaches  the 
"ablation”  temperature  at  about  200  kft. 

The  temperature  distribution  on  the  RV  surface  at  altitudes  5  seconds  apart  is 
shown  in  Figure  5  for  an  emissivity  of  0.25.  The  contour  at  the  crossover  between  free- 
molecule  heating  and  laminar  convection  is  indicated  by  the  dotted  curve.  The  steepest 
slopes  of  the  contours  along  the  surface  Oust  before  the  shoulder  on  the  hemisphere)  are 
only  2  to  3  percent  of  the  internal  temperature  gradients  just  within  the  surface  at  the 
shoulder,  shown  in  Figure  6;  two-dimensional  calculations  should  show  little  change  from 
the  one-dimensional  results. 

The  dependence  of  surface  temperature  on  emissivity  for  the  baseline  RV  and 
trajectory  at  four  altitudes  is  given  by  the  following  table  (in  degrees  Kelvin): 


emissivity  = 

0.25 

0.50 

0.75 

1.00 

stagnation-point  temp  at  295.8  kft  (90.17  km) 

886 

876 

866 

857 

stagnation-point  temp  at  266.3  kft  (81.18  km) 

1857 

1725 

1626 

1542 

conical-surface  temp  at  246.6  kft  (75.18  km) 

1021 

1001 

984 

969 

conical-surface  temp  at  207. 1  kft  (63. 13  km) 

1890 

1746 

1654 

1592 

The  fairly  weak  dependence  of  surface  temperature  on  emissivity  indicates  that  reradiation 
plays  a  smaller  part  than  conduction  in  dispersing  the  aerodynamic  heating  from  the  RV 
surface,  even  with  the  low  conductivity  of  glass-fiber  phenolic. 
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Flour#  4.  Temperature  History  on  ICBM  RV  Surface  at  Two  Points 
(aerodynamic  heating;  one-dlmenslonal  heat  conduction  In 
glass-fiber  phenolic;  emlsslvlty  «  0.25)  p  -  1500  lb/ft2 
(nose  radius  ■  0.077m);  y  *  -24.8°;  Rimpaet  *  10,020  km 
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Figure  5.  Temperature  Distribution  on  ICBM  RV  Surfacs  at  Dlffarent  Altitudes 
(asrodynsmle  heating;  one-dimenslonal  heat  conduction  In 
glass-fiber  phenolic;  emlssivlty  •  0.25) 
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Figure  6.  Temperature  Profile  Within  Heat  Shield  at  Hemisphere* Cone 
Shoulder  at  Different  Altitudes  for  the  Example  Reentry  Trajectory 
for  the  Example  RV  Configuration  With  a  Glass-Fiber 
Phenolic  Heat  Shield;  emlsslvlty  «  0.25 


V.  OBSERVATIONS 


Inclusion  of  the  effects  of  transient  heat  conduction  produces  a  dependence  of  total 
emittance  on  emissivity;  total  emittance  would  be  independent  of  emissivity  if  conduction 
were  ignored. 

The  values  for  the  total  gray-body  emittance  (ecT4)  for  the  temperatures  in  the  table 
in  the  EXAMPLE  RESULTS  are 


emissivity 

Total  Emittance  (watt/cm2)  at 

=  0.25 

0.50 

0.75 

1.00 

stagnation  point  at  295.8  kft 

0.87 

1.67 

2.39 

3.06 

stagnation  point  at  266.3  kft 

16.86 

25.11 

29.74 

32.07 

conical  surface  at  246.6  kft 

1.54 

2.85 

3.99 

5.00 

conical  surface  at  207.1  kft 

18.09 

26.36 

31.84 

36.43 

As  emissivity  is  reduced,  the  total  emittance  (the  integral  over  the  gray-body 
spectrum)  decreases,  even  though  the  temperature  increases.  (If  conduction  had  been 
ignored,  ectT4  would  have  been  independent  of  e  and  would  have  been  just  oTeq4, 
i.e.,  equal  to  the  aerodynamic  heat-transfer  rate.)  The  energy  received  by  a  broad-band 
detector,  such  as  a  bolometer,  would  decrease  by  a  factor  of  two  or  more  as  the  emissivity 
is  reduced  from  1.0  to  0.25.  For  a  detector  sensing  in  a  wavelength  band  beyond  the 
wavelength  of  greatest  spectral  emittance,  Xmax  (given  by  Wien's  displacement  law: 
^max  T  »  0.3  cm  K;  e.g.,  Xmax  »  3  p.m  for  T  =  1000  K),  the  ratio  of  the  energy  in  that 
band  to  the  total  energy  would  decrease  as  the  temperature  is  increased,  so  the  spectral 
emittance  (watt/(cm2K))  would  decrease  even  faster  than  the  total  emittance  as  the 
emissivity  is  reduced. 

On  the  basis  of  this  argument,  an  RV  designer  desiring  to  minimize  the  detectability 
of  his  RV  would  choose  a  surface  coating  with  the  lowest  practicable  emissivity  up  to  the 
temperature  at  which  the  surface  starts  to  break  down  and  evolve  gases. 
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APPENDIX 


TRANSIENT  HEAT  CONDUCTION 


A*1 


IMPLICIT  DIFFERENCING  TECHNIQUE 
(adapted  from  Richtmyer  and  Morton,  1967) 


Equivalent  Difference  Equation 


so  the  new  temperatures,  at  the  next  (n  +  1st)  time  step  in  terms  of  values  at  the  nth  time 
step,  are 

(3)  a  1^11-(2a+l)T^1  +  a  . 

If  we  assume  there  exists  a  recursion  formula  for  Tj+i  such  as 

(4)  T*+1  =  T”*1  +  d"  (compute  T  in  ascending  j) 


A-2 


then  (3)  becomes 


-  (2a+l)  T^+1 


c" 

j 


or 


so  the  recursion  formulas  for  c,  d  are 


(compute  c,  d  in  descending  j) 


Boundary  Conditions 

back:  Tu  =  cn  T10  +  du  =  Tl0  (insulated  wall  at  j  =  10|-) 


-*  cio  ~  a/(a+l)  ;  d®0  =  -^0  /a 


front 


T  _  HAx 
k 


^n+l  _  C^HAx/k  -d" 
1  = 


Tw  =  (3TQ  +  6Tj  -  T2)  /8  (quadratic  interpolation) 


Verification 

The  implicit  differencing  numerical  procedure  used  in  the  TREDE  program  (copy 
appended)  was  checked  against  an  analytical  solution*  to  the  diffusion  equation,  for  a 
constant  heating  rate  H,  a  thickness  b,  and  a  diffusivity  D  =  k/cp,  which  is 


0 


Supplied  by  Dr.  Irvin  W.  Kay. 


A-3 


T(x,  t)  =  T0  +  |-  Dt/b  -  b/6  +  (x-b)2/(2b)  -  ^  £(-!)" 

n  n=l 


nJt(x-b) 


The  procedure  with  a  Ax  =  b/10  (as  above)  and  a  At  =  0.1  sec  reproduced  the  values  given 
by  the  analytical  solution  to  within  1  degree  (out  of  <  300)  at  all  points  at  all  times.  For  a 
At  =  1  sec,  it  became  obvious  that  the  temperatures  given  by  the  numerical  procedure  were 
0.5  sec  behind  the  analytical  solution  (with  the  same  ~  1  deg  accuracy);  the  initial  time  step 
is  At/2.  Changing  the  Ax  to  b/100  from  Ax  =  b/10  for  At  =  1  sec  did  not  change  the  wall 
temperatures  by  more  than  1  degree  after  the  first  time  step. 


A-4 


I 


PROGRAM  TRIOE  ! RV-HEATING  W/WALL  TEMP  CORR  8/1 S/90 

DIMENSION  T(12) ,  C(12),  0(12),  P<12> 

DIMENSION  TEQ(44),  ALT (44) ,  VEL(44) 

DATA  TEQ/130., 142. ,134., 168. ,183., 200., 218., 23?., 239., 282.,  308., 
&  333.. 366.. 399. ,433., 473. ,323. ,376. ,636., 711.,  793. ,889. ,999. , 

&  1127., 1264., 1414., 1373., 1729., 1883.. 2039.,  2186. ,2324. ,2439. , 

&  2392., 2724., 2836., 2992., 3133., 3281., 3436.,  3601 ., 3788. , 3993 . , 

&  4214./ 

DATA  hLT/547. 5, 337. 9, 528. 4, 318. 8, 309. 2, 499. 6, 490. 0,480. 4, 470. 8, 

&  461.1,  451.5,441.8,432.2,422.3,412.8,403.1,393.4,383.7,374.0, 

&  364.3,  334.3,344.8,333.0,323.2,315.4,303.6,293.8,286.0,276.2, 

&  266.3,  256.3.246.6,236.8,226.9,217.0,207.1,197.2,187.3,177.4, 

&  167.3,  137.3,147.6,137.7,127.7/ 

DATA  UEL/22960 . , 22973 . , 22986 . , 22998 . , 23011 . , 23024 . , 23037 ., 23030 . , 
&  230  62 . , 230  75 . ,  23088 . , 231 0 1 . , 231 1 4 . , 231 27 . , 231 40 . , 23133 . , 231 66 . , 
&  23179. ,23192. ,23203. ,  23218. ,23231 . , 23244 . ,23237 . ,23270 . ,23283. , 
&  23296 . . 23310 . , 23322 ., 23336 . ,  23348 . , 23361 ., 23374 . , 23386 . , 23398 . , 
&  23409. ,23419. ,23428. ,23433. ,23439. ,  23440 . ,23434 . ,23419. ,23388./ 

102  FORMAT! 4X,  3HALT,  5X,  2HTW,  4X,  2HT1 ,  4X,  2HT2,  4X,  2HT3,  4X, 

&  2HT 4 ,  4X,  2HT3 ,  4X,  2HT6,  4X,  2HT7,  4X,  2KT8,  4X,  2HT9,  4X, 

&  3HT10 ) 

PI  -  3.1413926536 

SI G  *  3.672E-12/4.186 

EM  -  0.23 

XK  «  0.000826 

CP  «  0.3 

RHO  »  1.825 

TA  «  2000. 

TO  «  300. 

X  «  0.3 
OT  »  1. 

DX  -  X/10. 

103  FORMAT ( 2X ,  11H  INPUT  HFAC) 

1  CONTINUE 

TYPE  103 
READ! 5  ,*)  HFAC 
I F ( HFAC  .EQ.  0.)  GO  TO  99 
TYPE  102 
DO  I  >1.  12 
T( I )  -  TO 
END  DO 
TW  »  TO 
CC12)  *  1. 

0(12)  »  0. 

DO  10  K  -  1,  44 
TWX  «  TW 

IF(  K  .EQ.  3)  TWX  «  2.  *  TU  -  TUP 

I F( K  .EQ.  4)  TWX  *  3.  *  TW  -  3.  *  TWP  +  TWPP 

I F( K  .GT.  4)  TWX  ■  4.  *  TW  -  6.  *  TWP  +  4.  *  TWPP  -  TWPPP 

I F (TWX  .GE.  (TA-30 . ) )  GO  TO  1 

TWC  -  TWX 

Q  »  SIG  *  HFAC  *  (TE0(K))**4  * 

&  (VEL(K)*VEL(K)  -  6.73E6  *  TWC/289 .  )/(VEL<  K)*VEL(  K)  -  6.73E6) 

H  -  (-0  +  EM  *  SIG  *  TWX** 4)  *  OX/XK 
DO  I  -  1,  10 
J  »  12  -  I 

A  ■  XK  *  DT/(CP  *  RHO  *  2.  *  DX  *  DX) 

P(J)  -  -(T(J)  +  A  *  (T( J+l )  -  2.  *  T(J)  +  T(J-l))) 

C(J)  ■  A/(2.  *  A  +  1.  -  A  *  C< J+i ) ) 

D(J)  »  C(J)  *  (A  *  D( J+l )  -  P(J))/A 
END  DO 

T( 2)  «  (C(2)  *  H  -  D( 2) )/( C( 2)  -  1.) 

T( 1)  «  (T(2)  -  D( 2) )/C( 2) 

DO  J  »  3,  12 

T(J)  ■  C(J)  *  T(J-l)  +  D(J) 

END  DO 

TWPPP  -  TWPP 
TWPP  »  TWP 
TWP  «  TW 

TW  ■  (3.  *  T<1)  +  6.  *  T( 2)  -  T( 3) )/8. 

100  FORMAT ( 2X ,  F6.1,  2X,  11(F6.0)> 

TYPE  100.  ALT ( K ) .  TW.  (T(M) ,  M  ■  2.  11) 

10  CONTINUE 
GO  TO  1 
99  STOP 
ENO 


t 


